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Abstract 

As a model for a deformed nucleus the many level pairing model (picket fence 
model with 100 levels) is considered in four approximations and compared 
to the exact solution given by Richardson long time ago. It is found that, 
as usual, the number projected BCS method improves over standard BCS 
but that it is much less accurate than the more sophisticated many-body- 
approaches which are Coupled Cluster Theory ( CCT ) in its SUB2 version 
or Self-Consistent Random Phase Approximation (SCRPA). 
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I. INTRODUCTION 



The importance of two nucleoli pair correlations in the ground state and low lying excited 
states of nuclei has been known for a long time Q. The application to nuclear systems |2| 
of the concepts used in the description of superconductivity in solids was made immediately 
after the BCS theory has appeared ||. During the sixties it was realized that the pair- 
ing interaction was relevant in the description of two particle transfer reactions in normal 
and superconducting nuclei @J5|]. This interaction, which is usually thought to represent 
the short range part of the bare nucleon interaction, was treated by many authors in a 
phenomenological and schematic way. Nevertheless, it has been found recently |J that the 
pairing interaction is an important ingredient of the shell model interaction derived from 
realistic forces and used in large scale shell model calculations (the other two important 
ingredients being the quadrupole-quadrupole and the monopole-monopole interactions). It 
is also known J7| that to preserve the short range character of the force, it is necessary to 
use a large number of shells. Unfortunately the most simple theory for pairing in finite nu- 
clei, namely the mean field BCS approach, is rather limited in its application, since particle 
number fluctuations are very strong. Therefore more sophisticated approaches such as par- 
ticles number projection or the explicit introduction of quantal fluctuations like BCS-QRPA 
approach || and other more elaborated theories have to be considered. 

In a series of papers between 1963 and 1968 Richardson |§ obtained the exact solution 
of the pairing hamiltonian providing an analytic form for the eigenvalues and eigenvectors. 
These papers have found a revival in the framework of ultrasmall metallic grains [|H]] where 
it was necessary to go beyond the existing approximations to explain the disappearance 
of superconductivity as the size of the grain [0,0 decreases. Subsequently, the exact 
solutions have been generalized |13| and applied to other systems like Bose condensates |H 



interacting boson models ]15[ and nuclear superconductivity IfLql . 

The purpose of the present paper is to test on the exact solution for a large scale case 
the precision of some well known approximations like number projected BCS (PBCS) |n 



Coupled Cluster Theory (CCT) [H,and Self-Consistent RPA (SCRPA) @. In reality 
the possibility of applying these approximations depends on details of the nuclear residual 
interaction. In general these approximations can deal in an appropriated way with the long 
range part of the interaction, that can be thought of, in a simplified way, as a particle-hole 
interaction (as for example the quadrupole-quadrupole one), but they may have problems 
in dealing with the short range part that can be represented by the pairing interaction. 
The possibility or convenience of using each of these methods depends on the strength of 
the pairing interaction as well as the set of single particle levels that is considered. For 
example for very strong pairing (which is equivalent to a single shell) it is known that all the 
particles participate in the ground state wave function, and therefore one will need a quite 
large number of particle-holes over the Hartree-Fock (HF) groundstate to describe properly 
the paired state, on the other hand for a weak pairing interaction the ground state wave 
function will be almost the HF one. 

The paper is organized as follows. In Sect. II we present the picket fence model and 
sketch the main steps for its exact solution as given by Richardson. In Sect. Ill we outline 
the various approximate methods to treat the model and in Sect. IV we give the results 
together with a discussion. We end with the conclusions in Sect V. 



II. THE MODEL 



The picket fence model mimics superfluid correlations in a deformed nucleus where the 
level density can be considered as more or less constant and the levels are two-fold degenerate 
(for one sort of nucleons). As mentioned in the introduction, the model has been solved 
exactly in the early sixties by Richardson |§ for practically any number of levels. The latter 
feature makes the model very interesting because one can treat situations, very frequent in 
practice, which can not be mastered by ordinary diagonalization techniques. For example 
we here will treat the case of hundred particles distributed in hundred levels corresponding 



to a dimension of the hamiltonian matrix of 10 29 , well beyond any diagonalization technique. 
The model has not been used very much in nuclear physics, probably because of its schematic 
character. However, recently, its properties have been exploited in rather great detail in the 
context of ultra small superconducting metallic grains [p]Q[ . We here will employ the model 



in order to assess the quality of commonly used approximation schemes for nuclear pair 
correlations. We will consider fermion creation a) am and annihilation a am operators defined 
in a discrete basis labelled by the quantum numbers {am}. This basis can be referred to 
the single particle states of an external potential, the single particle energies depend on the 
quantum numbers a,m. 
The three operators 

n am — a a m a am i ^am ~ a am a oim ~ {A-am) (1) 

close the commutator algebra 



n ami Ap n 



2,So,fjS mn A.^ am , A am , Ag n — S a ^5 mn (1 n am ) (2) 



In Eq. (|l]) the pair operator A^ am creates a pair of particles in the time reversal states 
{am, am} where a/ am creates a particle in the time reversed state of a\ m . We will work 
with nucleons interacting via a pure pairing force and for simplicity we will represent by a 
single letter k the quantum numbers am (and when there is no possibility of confusion it 
will represent the pair {am, am}). Therefore the Hamiltonian that we will consider is 

H = Y,e k n k + GY,A\A kl (3) 

k kk' 

where the Sk are the single particle energies. 

The exact solution of this model has been obtained long ago by Richardson || . We will 
here briefly outline the method, giving the equations to be used later on in the numerical 
applications. 

Richardson || has shown that the exact eigenstates of the Hamiltonian (|3p with M pairs 
can be written as 



M 

i*> = rKi^> w 

i=i 

where there are v unpaired nucleons. The state \ip u ) describing the unpaired sector of 
is defined by the action of the operators A and n as 

A k \<p v ) = , n k \<p v ) = v k \ip v ) (5) 

where v k = 1 if there is one particle blocking the state k and v k = elsewhere. 
The operator B\ in (d) creates a collective pair 

where f2 is the number of single particle levels in the valence space. The form of the 
amplitudes in (|6|) were suggested by the one pair diagonalization of the pairing Hamiltonian 
(0). The pair energies Ei are unknown parameters to be determined by the eigenvalue 
condition. 

H\V) = E\V) (7) 

After a long but straightforward derivation one arrives at the set of M nonlinear equations 
for the M pair energies 

1-2G V -^ + GV (1 + 2i/t) =0 (8) 

j&i)=l ^3 k=l Zfcfc ^ 

while the energy eigenvalue is 

m n 

E = J2Ei + J2 £ ^k (9) 
i=i k=i 

The pair energies Ei are the roots of the set of M coupled equations There are as 
many independent solutions as states in the Hilbert space of M pairs. The different solutions, 
each one corresponding to an eigenstate of the pairing hamiltonian, can be classified in the 
limit of G — > as the different possible configurations of M pairs in Q levels, and then let 
them evolve adiabatically by solving the equations (H) for increasing values of G . 



The occupation probabilities are obtained by means of the Gellman-Feynman theorem, 
minimizing the energy with respect to the single particle energies 

dE 

ni = arr Vk+ ^w k (10) 

Differentiating (g) with respect to the occupation numbers can be expressed as 

i=i l^fc - A 

where the .Dj should satisfy the system of equations 

1 



2^' 



(11) 



^ (l + 2u k ) 

fc=l ( 2£ fc - ^i) 
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1 



(12) 



The above equations are used to establish the exact solution with, in the case considered 
here, a hundred levels with a hundred of particles. 



III. APPROXIMATE SOLUTIONS 



We will study some approximations that are written in terms of particular particle-hole 
excitations on a reference HF state. For simplicity we will consider the case when the shells 
are half filled, i.e. the number of pairs of particles M will satisfy Q = 2M. In the weak 
interaction limit the separation between the energy levels is much greater than the gap. The 
physics of this regime can be given in terms of the fluctuations around the HF state 

M 

i^> = n^io> (is) 

h=\ 

where h (p) refers to single particle states that are occupied (unoccupied) in the limit 
G = 0. 



A. Variational treatments 



We will first consider different variational treatments. The simplest one is the standard 
BCS treatment. The next approximation that we will consider is the number projected 



(before variation) PBCS wave function where the ground state is assumed to be a condensate 
of pairs of fermions. It is written as 

1 r_,lM 



where 



\PBCS) 



10} 



n m a 

r+ = £A^+ = £A^+ + £ \ p A; = r+ + r; 

k=l h=l p=M+l 



(14) 



(15) 
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(16) 



and |0) is the vacuum for the creation operator of the nucleons. In general one minimizes 
the energy by changing the variational parameters A&. In PBCS A& can be written in terms 
of the v k and Uk parameters as = with v l + ul = 1 . 

In Ref. |12| the ground state energy was evaluated in terms of these \ k coefficients using 
the auxiliary quantities 



j N,n 
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(17) 
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(18) 
(19) 
(20) 



and 
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The ground state energy is then written as 



(21) 



E gs = 2M £(2e, - /i) A.Sf + G £ A.Sf - GM{M - 1) £ \% 



2rpM 

ij 



(22) 



Kl 



The auxiliary coefficients are determined by recurrence relations using the fact that 
Z = l;Z 1 = Z l \l andSf = 

The pair creation operator has two parts: one (r+ ) creates two particles above the Fermi 
sea while the other part (T^ ) creates two particles below the Fermi sea. In Ref. [llj it is 
shown that if one defines the normalized states 

\K) = y^—(T;T h ) K \HF) (23) 

it is possible to write down the PBCS state as 

\PBCS >= ]T ^ BCS \K > (24) 

K 



where 



^,pbcs _ ((^/2)!) 2 Z K ,n/2 _ j^ Z k ,q/2 ^ 



and therefore the wave function can be written as 



PBCS >= An E V \HF) (26) 



(K\) 



For details on this derivation see fTT|| . 



The variational parameters in this wave function are the amplitudes It must be taken 
into account that An as well as the operators r+ and are well defined functions of these 
parameters. 

We also used another variational wave function with a structure similar to the exp{S-i) 
type(see below), i.e. 

K 



(r+r\) J 



\Exp >= B n £ V Vl \ HF ^ = Ba exp \ T p Th ) \ HF ^ ( 2? ) 

K A ' 

In this case the dependence on the parameters appears through the structure of T+ 
and and also in an indirect way in the normalization constant Bn- 



B. The Coupled Cluster Theory 



The CCT has been proven in the past to be a highly performant method for the cal- 
culation of correlation functions [18|. It has, however, never been tested for pairing model 
hamiltonians which is an interesting study case because of its exact solvability, even for very 
large number of particles. 

The Hamiltonian of the picket fence model can be written in the particle-hole basis as 

H = E &p + E frh -g|e 4 A p> + E 4 A' + E [ A l A h + A 1 A p] } (28) 

P h (P^P' h ^ h ' P h ) 

where 

e' p = e p — G/2 and e h = e h — G/2 



The unnormalized CCT wave function in the SUB2 approximation [M^ is 



ph 



(29) 



where the HF Slater determinant is given in (p~3|) . We have stopped at the one p-pair one 
h-pair, i.e. at the SUB2 level for reasons given below. The aim of the CCT is to determine 
the parameters x p h and the ground state energy. Acting with the Hamiltonian on the wave 
function we have 

H\m) = E\m) = Ee S2 \HF) (30) 
The key point of the CCT is to multiply (^) with e~ S2 from the left. Then 

e~ S2 H\m) = E\HF) (31) 



Projecting on the HF bra 



E = (HF\ e~ S2 He S2 \HF) , 



(32) 



taking into account that 



Si \HF) = (HF\ S 2 = 



27) is reduced to 



E = (HF\ He Si \HF) 
Having in mind the form of the Hamiltonian (|2q) , the groundstate energy is 

E = E HF — G Y x ph 

ph 

The amplitudes x p h are determined from the set of equations 

(HF\ A[A p e- S2 He S2 \HF) = 

which follows immediately after (^TJ) 
Eq. ( j36|) can be expanded as 

(HF\ A\A V (1 - S 2 ) H(l + S 2 + Sl/2) \HF) = 

The different terms are 

(HF\ A[A p H \HF) = -G 



(HF\ A{A P (-S 2 ) H \HF) = -x ph E HF 



(HF\ A{ x A{ 2 A hs A hi \HF) = (1 - 6 hlha ) (6 hlh4 5 h2h3 + 6 hlh3 8 h2h4 ) 



\HF\ A[A P HS 2 \HF) = 2 (e' p - e' h ) x ph + 2x ph Y e' h , - G 



Y X P'h + Yl X P h ' 



{HF\ A\A P (-S 2 ) HS 2 \HF) = Gx ph Y x P 'h' 



l -(HF\A{A p HS 2 2 \HF) = -G i 



p'h' 



x ph ^ t x p'h' ^ ] x ph' x p'h 



And therefore it is possible to write the equation for x p h as 

2 (e' p - 4) x ph + 2Gx ph Y x ph' + 2Gx ph Y x v'h + 2Gx 2 ph - G 



-GY X p'h ~GY X ph' + Gx ph -GY x p'h x ph' = 
p' h> p'h' 

This equation can be solved numerically. 



C. Self Consistent RPA 



The SCRPA for the Picket-Fence model has been developped in great detail in Ref. 

HH. Here 

we give a brief summary.The basic ingredients of the SCRPA approach in the 
particle-particle channel are the two particle addition operator 



h ■ 



(39) 



and the removal operator 



R{ = -Y, Y p X Qp + Y, x hQl 



(40) 



where Q p = A p /yl — (n p ) and Q h = —A' h /y(n h ) — 1. Where the expectation values are 
referred to the SCRPA vacuum defined as 



\SCRPA) = R x | SCRPA) = 



(41) 



and the collective RPA excitations are 
\N + 2)^ = A\\SCRPA) 



\N - 2) A = R[ \SCRPA) 



(42) 



The equation of motion method applied to these operators leads directly to the SCRPA 
equations 



(43) 
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where 



A pp/ = (0\[Q p ,[H,Ql,}}\0) 



5 pp i < 2e p + G + 2 
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-G 



1 - (np) 77 
((1 -n p )(l -n p /)) 

(i-K))(i-(v)) ' 



B ph = (0\[Q p , [H,Qi}}\0) = G j (1 1}) (44) 

^J(l - (n p ))((n h ) - I) 



C w = {0\\§ h ,[H,Ql,}}\0) 
= 8 h h 



+G- 



^/«n h >-l)«n fc ,>-l) ' 
Since the amplitudes X and Y form a complete orthonormal set of eigenvectors in (fl3| ) one 
can invert the Bogoliubov transformation of fermion pair operators ( |3l3[|4"0"|) and all expec- 
tation values in (fPj) can be expressed in terms of the RPA amplitudes and of the number 
operators expectation values (n p ), (n h ), (n p n p /), (n p n h ), and (n h n h i). For the particular case 
of the Picket Fence models these expectation values can be calculated exactly within the 
SCRPA approximation as shown in []21j ]. In this way the SCRPA constitutes a closed set of 
equations without any further approximation than the definition of the collective operators 
(^,^3) and the corresponding vacuum condition ([41]). 

Knowing these expectation values we can evaluate the SCRPA ground state energy: 

(H) = E4W + - G ( £ (AlA p ,) + £ (A[A h ,) + XX44i + 4A,) 1 (45) 

Assuming that the single particle energies are all equally spaced, separated by an 
energy gap e, we have e\ = Si — e/2 + G/2. In this case the SCRPA correlation energy is 

E^ PA = (H)+eM 2 . (46) 



D. Results and discussion 

We will study the approximate descriptions of the pairing interaction in the deformed 
nuclear region characterized by a constant density of levels near the Fermi surface. This 
situation, therefore, can be represented by a set of equally spaced levels with the appropriate 
density. We have used 100 levels with a constant level spacing of 300 keV and with 100 
nucleons (half filling). This represents typical values of the level density and neutron numbers 



in the rare earth region ( A ~ 170). As in this region the gap has a value of the order of 
A ~ 0.8 MeV the physical value of the pairing interaction G ~ 0.1 MeV. For this level density 
and number of particles the critical pairing strength of the model in the BCS approximation 
turns out to be G c ~ 0.055 MeV. 

The aim here is to compare the quality of different approximations to treat the pairing 
problem which are outlined in the text. We display in Fig. 1 the ground state energy 
obtained using the various methods discussed in the previous section (only the correlation 
energy is displayed to isolate the effects due to the interaction) . All the correlations energies 
are given in terms of the exact energy. Standard BCS approximation provides a rather poor 
description. The numerical results do not appear in Fig. 1 because they are out of scale. A 
strong improvement over BCS is obtained with the number projection before variation, i.e. 
the PBCS procedure. Still quite a bit better works the Exp method for moderate values of G 
, described at the end of section III. A, with the factorisable ansatz in the exponential. Both 
curves show a typical structure: for small G there is a linear regime which can be qualified as 
the perturbative regime. It is followed by a part with negative curvature, characterized by 
precritical fluctuations, before the superfluid regime develops after the minimum. A detailed 
study of the two former regimes has been performed in ref. p0[ . The figure also shows a 
clear indication that the PBCS approximation approaches the exact groundstate energy 
in the large G limit while this is not the case for the Exp method. Both approximations 
underbind, as it should be for a strictly variational theory in the sense of Raleigh-Ritz. 
On the contrary CCT (expS^) and SCRPA overbind because neither CCT nor SCRPA in 
general correspond to a Raleigh-Ritz theory. However, in absolute values both of the latter 
theories work extremely well. It should be pointed out that since SCRPA is a theory for 
two body correlation functions, we only can go in CCT up to the SUB2 approximation, for 
consistency. Going to higher approximations, we should also include higher than two-body 
correlations and SCRPA and CCT would not be on the same level of approximation. 

We only have worked in the normal particle basis for CCT and SCRPA and therefore the 



iterative solution of the eqs (^) and (|44|) did not converge any longer beyond G/G c ~ 1.3. 



We know from experience in other models |T9| that around the mean field phase transition 
point one has to change to the "deformed" basis which means to the quasiparticle basis in 
our case. For PBCS and Exp the error in the correlation energy in the superfluid phase 
decreases for G > G c and therefore the correlation energy has its maximal error in the 
transition region as it is to be expected. For the picket fence model we have not yet worked 
out the SCRPA in the superfluid phase and we are not aware of any attempt to apply CCT 
in this regime. As mentioned before, both curves in Fig. 1 stop at the point where we 
do not find a numerical solutions of the corresponding equations any more. We, however, 
conjecture that the end points of both curves represent the maximal error and continuing 
the calculation in the superfluid phase the error would start decreasing again. We see that 
the errors in expS2 and SCRPA are, in the worst case, only of 5% and 2% respectively. 
These errors are much smaller than PBCS and Exp which are of the order of 15% — 20%. 
The very small errors of expS2 and SCRPA is a very satisfying result which confirms earlier 
positive results with these theories for correlation functions in other cases. The factor two 
improvement of SCRPA over expS2 for the correlation energy in the transitional region 
has already been found in another model study JL9J but this may be accidental. Grossly 
speaking both methods are of similar characteristics and accuracy for the correlation energy 
in the normal phase. The main advantage we see in SCRPA is that excitation energies 
and correlation functions are obtained simultaneously from the same theory. The SCRPA 
excitation energies also turn out to be very accurate in the present model (see ref. [[HJ]). In 
CCT the excitation energies have to be constructed separately putting new ingredients into 
the theory. 

In conclusion in this work we have compared four methods for the calculation of energies 
in the pairing case with parameters typical for deformed nuclei. This study was performed 
in the picket fence model with a model space of a hundred levels. The exact solution could 
be obtained owing to the method proposed by Richardson long time ago, whereas a brute 
force diagonalization is far beyond the limits of present computers. We found that the expSi 
and the SCRPA methods are quite superior to the other variational methods in the normal 



phase. The results obtained in this work might stimulate further efforts to extend both 
approximations to the superfluid regime and more realistic forces. 
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Figure Captions 



Figure 1: Ratio between the approximate and the exact correlation energies for equally 
spaced levels as a function of the pairing strength for the four approximations discussed 
in Section III. 
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